Principles of Numerical weather prediction and modeling 


There are three types of weather forecasts. They are short range, medium range and long range 
weather forecasting. 

Short range = few hours to 48 hours or sometimes 72 hours 

Medium range = 3 days to about 3 weeks 

Long range = a few seasons to a year 

There are three weather forecasting methods. They are i) synoptic weather forecasting ii) Numerical 
methods and iii) statistical methods. 

Synoptic forecasting: 

Synoptic means the observation of different weather elements refers to a specific time of observation. 
From the careful study of weather charts over many years, certain empirical rules can be formulated. 
These rules will help to estimate the rate and direction of movement of weather systems like a trough 
line, hurricane or thunderstorm. 

Numerical methods: 

Numerical method is based on the fact that gases of the atmosphere follow a number of 
physical principles. If the current conditions are known, these physical laws help to forecast the future 
condition. 

A series of mathematical equations are used to develop theoretical models of the general 
circulation of the atmosphere. These equations are used to specify changes in the atmosphere as time 
changes using the weather elements like wind, temperature, humidity, clouds, rain, snow, evaporation 
etc. 

The atmosphere is divided into 6 or 11 layers for this purpose. As the initial conditions ( initial 
state of the atmosphere) are known through observations using surface and upper air data, it is possible 
to predict using numerical models for future conditions. 

Valuable contribution in the development of numerical methods was made by Prof.J.Charney in 
1948 and Obukov. 

As there is large number of weather variables and only a less number of equations are available, 
some weather variables are ignored on the assumption that they do not change with time. As large 
numbers of calculations are involved even with simplified models, a high speed super computer is 
necessary. 

For carrying out numerical models the following considerations are followed: 

The equations are simplified 

The equations are designed in such a way that they guarantee the conservation of air mass, momentum, 
water vapor and total energy for all time. 

The region to be forecasted is divided into number of nodal points in a rectangular shape. 

The equations are solved at each nodal point for a short period of 10 minutes. 

By repetitive calculations for every next 10 minutes forecast is obtained for 24, 48 or 72 hours. 

L.F. Richardson in 1922 attempted to solve the differential equations of momentum, energy, state and 
conservation using finite difference method with the help of desk calculator. 

After World War II and the development of digital computer, it was realized that Richardson's 
scheme was not simple as the governing equations contains not only slow moving meteorologically 
important motions but also high speed sound and gravity waves. Though these waves are very weak, 
they amplify surprisingly in numerical methods and thus generate 'noise'. 

So in 1948 J.G.Charney simplified the equations by introducing the scale analysis, geostrophic 
and hydrostatic assumptions. By these techniques sound and gravity waves are filtered. Then these 


equations are called quasi-geostrophic models. A special case if this model , called equivalent barotropic 
model, was used in 1950 which gave the first successful numerical forecast. 

A. Objective Analysis 

The automatic process of transforming of irregularly distributed observed data from different 
geographic locations at various times to numerical values at regularly spaced grid points at fixed times is 
called objective analysis. In the earlier days of numerical weather prediction, the two-dimensional fitting 
of polynomial was applied to observational data in an area surrounding a grid point at which the value 
of analysis is required. Polynomial fitting is suitable for processing relatively dense and redundant data, 
but it often gives unreasonable results in data of sparse areas. 

More recent procedures of objective analysis practiced at most of the operational forecasting centers 
consist of the steps schematically described in Figure below. Most observed data are collected 
worldwide every 6 hours, say 00, 06,12 etc. UTC (Coordinated Universal Time), and are blended into 
forecast values appropriate to the map time by the process referred to as data assimilation. The basic 
process of data assimilation, indicated by the box enclosed by dashed lines in Figure below, consists of 
two steps: One is objective analysis and the other is called initialization. The step of objective analysis is 
carried out by statistical interpolation by taking into account all of the available observations plus other 
prior information. Short-term forecasts from the previous analysis cycle are used as prior information 
and are referred to as background fields. 

B. Concept of Initialization 

The solutions of primitive equation models correspond to two distinct types of motion. One type has 
low frequency. 

Its motion is quasi-geostrophic and meteorologically dominant. The other corresponds to high- 
frequency gravity-inertia modes. The amplitude of the latter type of motion is small in the atmosphere. 
Flence, it is important to ensure that the amplitudes of high-frequency motions are small initially and 
remain small during the time integration when the primitive equations are solved as an initial value 
problem. The process of adjusting the input data to ensure this dynamical balance is called initialization. 
Solution to the initialization problem was central to the successful transition in forecast practice from 
the use of 

quasi-geostrophic models to primitive equation systems during the 1970s. Since gravity-inertial motions 
are filtered 

out in quasi-geostrophic models, no special procedure was necessary and the objectively analyzed data 
could be used immediately as the input data to quasi-geostrophic models. Actually, there is an 
interesting history that led to solution of the initialization problem that is equivalent to the quest of 
understanding what the primitive equation forecast system really is. A promising break came with the 
advent of so-called nonlinear normal mode initialization (NNMI) during the latter part of the 1970s. 
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Numerical solution of equations 


The five prognostic equations used for forecasting are: 
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The first three equations are equation of motion and the other two are equations of continuity and law 
of conservation of energy for an adiabatic process. There are five dependent variables u,v.w,p and p 
and five equations. So they can be solved. 


Finite difference scheme to solve the equations 


Using the Taylor's equation we can write 
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Taking first two terms in the right hand side, we can write 
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The finite difference grid is as below: 
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For example, to get ^ at middle point 'O' we can write: 

dp _ P2— P4 
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Similar expressions can be written for all other five variables. Thus incorporating all these values at each 
grid point, future time value is known. 

Filtering of Sound and Gravity waves 

Sound waves are longitudinal waves which are generated because of alternating adiabatic 
compression and expansion of the medium. This arises because of pressure gradient. In the case of 
hydrostatic atmosphere, the vertical pressure gradient can not be influenced by adiabatic compression. 
Hence the vertically propagated sound waves are filtered under the hydrostatic approximation. 
However, horizontally propagating sound waves are eliminated by the assumption that at the lower 
boundary co = = 0. 

After the sound waves are filtered, the prediction equations, such as momentum equation, continuity 
equation and hydrostatic thermodynamic equation are written as: 
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Here er = and the lower boundary condition is co = 0. The above set of equations along with the 
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hydrostatic and equation of state are called primitive equations. They contain the gravity waves as 
solutions. 

The gravity waves are waves in which buoyancy acts as the restoring force on air parcels 
displaced vertically. A disturbance is propagated because of alternate horizontal convergence and 
divergence in individual columns of fluid. A divergent horizontal velocity field which changes in time is 
essential for the propagation of gravity waves. Neglecting the local rate of change of the horizontal 
divergence in computing the balance between mass and velocity fields is sufficient to filter out the time 
dependent gravity waves. As the prediction equations of 1 to 3 do not contain the local derivatives V. V 
explicitly, we use the vorticity and divergence equations instead of horizontal momentum equations. 

The vorticity and divergence equations are as below: 
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The divergence equation is 
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By operating V. ( ) on this equation, the divergence equation can be written as : 
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If the Left hand side of equation (5) is set to zero, all solutions corresponding to time dependent gravity 
waves can be eliminated. But the equation (5) becomes an intractable diagnostic relationship between 
,V and (p . Hence we do the scale analysis for systematic elimination of unimportant terms. Usually 
neglecting the second order terms is sufficient to filter the gravity waves. 

In order to exhibit the connection between scaling and filtering, it is convenient to partition the 

wind as 

V = V xp + V e (6) 

Where is the non divergent part and V e is the divergent part of the wind. 

If the velocity field is two dimensional, the non-divergent part can be expressed in terms of 
stream function as 

V^=KxVip (7) 

As V. Vjp = 0, V X V e = 0 and 

$ = K.VxV = V 2 ip (8) 

The horizontal divergence is small compared to the vorticity (by scale analysis) for synoptic scale 
systems in mid latitudes. That is V is quasi- nondivergent 



Thus to a first approximation we can replace V by everywhere in equations 4 and 5except in terms 
involving horizontal divergence. 

So the resulting vorticity and divergence equations after scaling are 


f = -VV(C + /)—/„ V . Ve (9) 

V 2 cp = — /qV. (KxV, p ) = f 0 V 2 xp (10) 


In the above equations gravity wave terms were filtered. From equation (10), we can interpret as 
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This means the geopotential field on a constant pressure chart is approximately proportional to the 
stream function field. Equation 10 also says that to a first approximation, the vorticity is the vorticity of 
geostrophic wind calculated using the constant mid latitude value of the coriolis parameter. 

The geostrophic vorticity equation and the hydrostatic thermodynamic equation can be written 


as: 
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By eliminating a> from these equations, the quasi-geostrophic potential vorticity equation as: 

(£ + V„. V)q = 0 (13) 

Where , = (V 2 * + /) + / 0 2 ±Qg) 

We have assumed that a is a function of pressure only. This relation states that the geostrophic 


potential vorticity is conserved following the non-divergent wind in pressure coordinates. 

In practice solving equation (13) is extremely difficult and so one has to go for numerical schemes. 
Models have been developed which simplify the numerical structure by referring to one, two or few 
layers. 

One Parameter models 


There are some simplest models with only one cluster level in the vertical. 

The barotropic model: 

In this model one assumes that there exists a level of non divergence in the atmosphere. Usually this 
level is assumed as 500mb level. 

Thus V. V = — — =0 at 500 mb level 
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This equation allows to compute the evolution of the fluid at the level of non divergence only. 

b) The equivalent barotropic model: 

Same type of variation of wind with height is taken in equivalent barotropic model by averaging in 
the vertical. But one level is only considered in the model . In this model the nondivergent wind is 
expressed as: 


V, p (. x,y,p ) = A(p)< Vp, (x, y) > 

The angle bracket denotes the vertical average and A(p) is a weight function. The vorticity equation 
for this model is written at the stairred level where 

A(p ) = < A 2 (p ) > 


<■ and xp are also expressed as, 


<■ = i4(p) < > and xp = A(p ) < xp > (15) 

To the extent that equation (15) is correct, the results from the solution of vorticity equation can be 
utilized to deduce the evolution of xp at other levels. 


Two parameter baroclinic model 


The one parameter models do not allow temperature advection. But the thermal advection processes 
are essential to the development of synoptic systems. The barotropic models are quite effective in 
predicting the evolution of mid tropospheric flow for periods upto a couple of days. For consistently 
reliable forecast one should predict the development of new systems also. Hence there is a need of 
baroclinic models. The simplest baroclinic model is a two layer model. The atmosphere is divided into 
four discrete layers and the disposition of variable is given in the following schematic diagram. 
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Here T is level and p is pressure surface. 


For a flat ground the lower boundary condition is <d 4 = 0. Applying the vorticity equation (11) at 
levels 1 and 3 and utilizing the finite difference technique, we can write: 
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The vorticity equations at the two levels are: 
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Writing the thermodynamic equation (12) at level (2) by utilizing the approximation 
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Now equations 16, 17 and 18constitute a close set of prediction equations. 


First eliminate o 2 between 16 and 17 by adding both equations to get 
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This equation governs the barotropic part of the flow. 

Next subtract ( 17) from (16 )and add the results to V times (18) to get 
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Here A = — — is anon dimensional parameter. Equation (19) states that the local rate of 
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change of vertically averaged vorticity is equal to the average of vorticity advections at the two 
levels. Equation (20) gives the local rate of change of thickness of 250-750 mb layer. 

The time derivatives between equations 16, 17 and 18 can be eliminated to get the omega equation 
as 
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Since equations 19, 20 and 21 are derived from the quasi-geostrophic system of equations, in the two 
parameter model, the vorticity changes are geostrophic and the temperature changes are hydrostatic. 
Due to these constraints, the omega field is to be adjusted at every instant. 

Primitive Equation models 

The filtered models discussed above, involves various approximations of the quasi-geostrophic system. 
So there are some limitations on the accuracy of forecasting. In low latitudes, the stream function ip can 
not be obtained from the geostrophic vorticity field and can be obtained from more complicated 
balance equation which is a nonlinear equation. It is difficult and time consuming to solve the balance 
equation at each time step. 

The following balance equation is derived from the divergence equation by putting initially 
|-(V.y) = 0 and V.V = 0 . 

ot 
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Therefore, it is worthwhile to investigate the possibility of using less restrictive modeling assumptions. 
The scale analysis indicated that the hydrostatic approximation is far more accurate than other filtering 
approximations. Hence we continue to assume that the motions are hydrostatic. The resulting system of 
equations is isobaric coordinates consists of three prognostic equations (x and y components of 
momentum equations and thermodynamic equation) and three diagnostic equations (the continuity, 
hydrostatic ad the equation of state). 

These six equations constitute a closed set with the dependent variables u,v,w,a, 6, and 0 
which are referred to as primitive equations. 

Let us write the basic primitive equations in isobaric coordinate system as: 
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The transformation from p to a coordinates can be derived in a similar manner of transformation 
from Z to P coordinates. 


Isobaric and sigma coordinates 



In isobaric coordinates we put w(p 0 ) = — pogw(Z 0 ). We assume that the height of the ground 
Z 0 is coincident with the pressure surface p 0 ( = lOOOmb). However, this is not strictly true, as pressure 
surface will not coincide with the ground even when the ground is level. So there is a need of a system 
of coordinate is used in the vertical in numerical modeling such that a = — where p s denotes the 

Ps 

pressure at the surface. In sigma coordinate system the lower boundary condition is applied at a = 1 . 
The vertical velocity in a coordinates is a - will always be zero at the ground as a = 1 is constant 
even in the presence of sloping terrain. Hence in a coordinates the lower boundary condition is a = 0 at 
a = 1. 

Using the definition of , one can write: V„ ( ) = ( ) — — V„ — ( ) 
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So the system equations in a coordinates can be written as: 
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The set of six equations in (23) contains seven dependent variables u, v, a , 8, T <p and p s . So another 
equation is necessary to solve these variables. This equation is surface pressure tendency equation. It 
can be obtained by integrating the continuity equation vertically by using the boundary conditions a = 0 
at a = 0 and 1. 

dv 1 

The surface pressure tendency equation is = ~ J 0 (.PsV)d a 

This equation states that the rate of increase of the surface pressure at a given point is equal to the 
mass convergence in a unit cross section of column above that point. 

Two layer primitive equation model (PE Model) 

In the primitive equation model the static stability parameter a is a function of space and time 

d6 

whereas in earlier models, a (p) is usually a specified parameter. It can be shown that if — is specified 

as a function of pressure alone, that is a = a(p ), it is not possible to formulate a PE model which 

d6 

satisfies the requirements of energy conservation. Hence it is necessary to compute — explicitly at each 

dp 

time step. Hence temperature must be predicted at a minimum of two levels rather than at one level as 
in the case of Z-level quasi-geostrophic model. 



The schematic diagram of Z-level PE model with variables is as below: 



The momentum equation and the thermodynamic equations are applied at levels 1 and 3. The 
vertical differencing of these equations requires knowledge of sigma vertical motion at level 2. ri 2 is 
obtained diagnostically through the use of continuity equation as follows: 



From the continuity equation in (28) we can write: 

^ + V.(p s F 1 )+ 2p s & 2 = 0 (24) 

^ + V.(p s K 3 )- 2p s 6, = 0 (25) 

Adding (24) and (25) we can write: 
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Subtracting (25) from (24) we can get: 
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In a similar manner <p 1 and 0 3 can be obtained diagnostically by vertical integration of hydrostatic 
relation in equation (23). 



In short the following steps may be followed for the forecast using PE model. 

1. Suitable finite difference analogues of the momentum and thermodynamic equations at levels 1 
and 3 and the surface pressure tendency equation are written. 

2. These prognostic equations are utilized to obtain tendencies of Vj, V 3 , 0 1 , d 3 and p s . 

3. Extrapolating the tendencies ahead using a suitable time differene scheme. 

4. The new variables of V t , V 3 , 9 ll 8 3 and p s are used to determine d 2 , <fii and <p 3 diagnostically. 

5. Repeat steps 2 to 4 until the desired forecast time is reached. 

It must be noted that the PE model contains the mechanism for both gravity wave and horizontal 
acoustic wave propagation. It is necessary that the time increment be kept small enough to satisfy the 

condition: ^ where C is horizontally propagating sound wave which is the fastest moving wave 

solution, d is horizontal resolution and At is time increment. Therefore the time steps in PE models are 
considerably small ( = 10 minutes). 

Initialization Scheme 

Initialization scheme involves the use of accurate initial data. One reason for the failure of Richardson is 
due to lack of accurate initial data. Though the data coverage is vastly improved since then actual 
measured wind and pressure data is still lacking as most of the data is derived from remote sensing 
satellite data. 

To introduce the initialization method, consider the relative magnitudes of the various terms in 
the momentum equation in isobaric coordinates: 

dV dV 
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For synoptic scale motions in mid latitudes, we have seen that the wind and pressure fields are in 
approximate geostrophic balance. Thus the acceleration following the motion is measured by the small 
difference between the two nearly equal terms Coriolis and geopotential (—fK x V and Vcp ). Although 
cp can be determined with good accuracy using observed data, observed winds are often 10 to 20% in 
error. Hence if wind is used as initial data, the Coriolis force works out to 10 to 20% in error. Since the 
acceleration is normally 10% of the Coriolis force, the acceleration comes to 100% in error (10x10). Such 
spurious acceleration leads to poor estimates of initial pressure and velocity tendencies. This further 
produces large amplitude gravity wave oscillations. 

To avoid this problem, let us derive the balanced equation as follows. In the divergence equation 
at t = 0, put — + V. V = 0 and V. V = 0 
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Since we have assumed that the initial velocity is non divergent, we can replace V by 0 by letting 


= K x V0 then we get 

V 2 [0 + y 2 (V0) 2 ] = V. [(/ + V 2 0)V0] (26) 

This equation is called balance equation. Equation (26) is linear in 0 and non linear in 0. In midlatitudes 
0 is known and henc e 0 is solved with the help of equation 26. The 0 has two solutions. One solution 
corresponds to the gradient wind balance and the other one is anomalous. One 0 is known, the initial 
velocity field can be computed. This velocity is in gradient wind balance with the geopotential field. 
The initial local accelerations are thus very small. Since the initial rate of change of divergence is zero, 
the large amplitude gravity wave oscillations are not excited by the initial conditions. 



